Serially-regulated biological networks fully realize a constrained set of functions 
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We show that biological networks with serial regulation (each node regulated by at most one other 
node) are constrained to direct functionality, in which the sign of the effect of an environmental input 
on a target species depends only on the direct path from the input to the target, even when there 
is a feedback loop allowing for multiple interaction pathways. Using a stochastic model for a set of 
small transcriptional regulatory networks that have been studied experimentally fl] , we further find 
that all networks can achieve all functions permitted by this constraint under reasonable settings of 
biochemical parameters. This underscores the functional versatility of the networks. 



A driving question in systems biology in recent years 
has been the extent to which the topology of a biologi- 
cal network determines or constrains its function. Early 
works have suggested that the function follows the topol- 
ogy [H m [31 H] , and this continues as a prevailing view 
even though later analyses (at least in a small corner of 
biology) have questioned the paradigm |5j [6] . It remains 
unknown if a small biochemical or regulatory network 
can perform multiple functions, and whether the func- 
tion set is limited by the network's topological structure. 
To this extent, in this paper, we develop a mathematical 
description of the functionality of a certain type of bio- 
logical network, and show that the answer to both ques- 
tions is "yes": the networks can perform many, but not 
all possible functions, and the set of attainable functions 
is constrained by the topology. We illustrate these results 
in the context of an experimentally realized system ^ . 

Following [1] and our earlier work 6 , we focus on the 
steady-state functionality of transcriptional regulatory 
networks. In this case, the input is the "chemical en- 
vironment," that is a binary vector of presence/ absence 
of small molecules that affect the regulation abilities of 
the transcription factors; and the output is the steady- 
state expression of a particular gene, hereafter called the 
reporter. Different functions of the network correspond 
then to different ways to map the small molecule concen- 
trations into the reporter expression. 

In our setup, the effect of introducing a small molecule 
Sj specific to a transcription factor Xj is to modify the 
affinity of Xj to its binding site. Equivalently one can 
think of Sj as modulating or renormalizing the transcrip- 
tion factor concentration Xj by some factor Sj, making 



the effective concentration Xj = Xj i-^j T^j)- ^ simple 
example of such a modulation function is 

Xj{Xj,Sj) = Xj/sj, (1) 

in which the presence of the small molecule reduces the 
effective concentration of transcription factor by the fac- 
tor Sj. 

The function of the circuit will depend on how the 
steady-state expression G* of the reporter gene G changes 
as the modulation factor sj is varied from some "off" 
value sJ to some "on" value s^: 

AG* ^ G*{s+)-G*isJ) ^ J_ dG^^^ .2) 
As J Asj Asj dsj 

where Asj — ~ . For example, if Xj — ^j/sj, then 
sJ = I, indicating that the small molecule is absent, and 
> 1 is the factor by which effective concentration is 
reduced when the small molecule is present. 

If the sign of dG* /dsj does not change for sj G [sJ , Sj'], 
then the sign of AG* is fixed. For networks with only 
serial regulation, i.e. each gene is regulated by at most 
one other gene, we will show that the sign of dG* /dsj 
is unique and in accord with the direct path from Sj 
to G, a property we term direct functionality. This con- 
strains the possible responses and hence the functionality 
of serial networks. Importantly, we will then show that 
all admissible functions indeed can be attained by all the 
networks we studied operating at different parameter val- 
ues. While throughout this work we focus on the setup 
pioneered experimentally by Guet et al. pj , we also show 
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FIG. 1: Four-gene networks (3 transcription factors Xi plus 
1 reporter gene G) in which each gene is regulated by one 
other gene, as studied in [T[|2]. Transcription factor efficacies 
are influenced by small molecules Si. Regulation functions Qi 
are assigned to the edges. The three edges ai, 32, and as 
can be up-regulating or down-regulating, giving 3 x 2^ = 24 
possibilities; the reporter gene is repressed in all cases. 



that the constraint to direct functionality holds for any 
network with serial regulation. 



DIRECT FUNCTIONALITY IN SMALL 
NETWORKS 



As in Guet et al. [T], we consider networks with N = 4 
genes (three transcription factors plus a reporter G), in 
which each gene is regulated by exactly one other gene. 
This admits three topologies and a total of 24 networks, 
as described in Fig. [l] All three topologies consist of a 
cycle and a cascade that begins in the cycle and ends at 
the reporter gene G. Once outside the cycle, there is only 
one path to G, so it suffices to study a topology consisting 
of an n-gene cycle with a gene G immediately outside 
(Fig. [ij; is an example with n = 3), and extensions to 
topologies where the cycle is connected to the reporter 
by a linear cascade are trivial. 

In this section, we will perform the steady-state anal- 
ysis of such single-cycle networks to lay the groundwork 
for understanding the effect of topology on allowed func- 
tionality. 

The process of protein expression has been modeled 
with remarkable success by combining transcription and 
translation into one step and directly coupling genes by 
a deterministic dynamics [H HI |^ . Accordingly we model 
mean expressions Xi (we later distinguish between entire 
probability distributions P{Xi) and the means of these 
distributions Xi] cf. Appendix) with the system of ordi- 



nary differential equations 
dXi 



dt 
dX, 
dt 
dG 



aiixn) -riXi, (3) 
a,{x^-l) - nX, {2<i<n), (4) 

an+l{Xn) - Tn+lG, (5) 



where the ai are creation rates for the species X.^ (and 
Xn+i = G), each monotonically regulated by the effec- 
tive concentration XTTi of its parent tt^, and the r; are the 
decay rates. Note that we have set 



^ i-1 



{2<i <n + l) 



(6) 
(7) 



to create the n-gene cycle with one gene immediately 
outside. The regulation functions cti will be up- or down- 
regulating according to the network topology. A common 
example is the familiar Hill functions. 
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(up-regulating) , (8) 
^ (down-regulating) , (9) 



with basal and maximal expression levels Uq and Oq + a 
respectively, Michaelis-Menten constants K, and cooper- 
ativities h. Although we use the functional forms in Eqns. 
([8]|9|, as well as the functional form for the modulation 
function in Eqn. ([T|), for our numerical experiment (cf. 
Numerical Results), the analytic result derived in this 
section will be valid for any monotonic functions a{x) 
and any function xi-^^ •s)- 

Fixed points of the dynamical system in Eqns. (|3][5]) 
satisfy 

X*i = ai(x:j, (10) 
X* = a,(x*_i) (2 < z < n), (11) 

G* - «„+i(x:), (12) 

where we define 

ai = ai/ri. (13) 

We may now, as in [101 Eli use the chain rule to cal- 
culate the derivative of G* with respect to a particular 
input factor Sj. For illustration, we will do so first for 
the concrete example in Fig. [Tj:, in which n = 3. Let us 
consider the derivative of G* with respect to si: 



dG* 
dsi 



da4 das 



da2 
dsi 



da2 dXl 



dXi dsi 



(14) 



where all derivatives are evaluated at the fixed point, 
and it is understood that ai depends on either X„. or 
STTi through x-Ki, that is, that 

dai dui dxTTi , d^i da^ dxi^ 
— and — 



dX^^ dx-Ki dX^^ 



(15) 
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If we introduce the notation 



then Eqn. (14 1 becomes 



dG* 
dsi 



, dxi 

"2 + Oi2—, 



(16) 
(17) 



(18) 



The first term reflects the direct chain to G from Si, 
and the second term incorporates further contributions 
around the cycle and will need to be evaluated self- 
consistently. 

For a cycle of arbitrary length n and for an arbitrary 



input factor Sj (1 < i < n), Eqn. (18) generalizes to 



dG* 
ds-j 



dx; 



n+1 

n " 

J k=j+2 



(19) 



where we use the convention that 

b 



Jl [•] = 1 if a>b. 



(20) 



k—a 



We may also use the chain rule for dX^ / dsj , 



d_xi 

dsi 



Q!(jmodn) + l + Q!(jmodn) + : 



dSj 



nn I 



(21) 



[j mod n) + l 

and now we may solve for dX* /dsj self consistently: 
dX* _ d(jmodn)+i nLi"fc 



(22) 



For the special case of j — n, where (j mod n) + 1 = 1 



substituting Eqn. (22 I into Eqn. (19 I obtains 



dG* 

dSn 



1 



n 

1 



fe=2 



where the second step follows from 

. , f dai dxn \ f dan+i dxn 



\dXn dSn J V '^Xn OX^ 

f dai dxn \ ( dan+\ dXn 



\dxn dX„ 



dXn dSn 



(23) 
(24) 

(25) 

(26) 
(27) 



in which the first step recalls Eqn. ( 15 1 . For 1 < j < n— 1, 



where (jmodn) + 1 = j + 1, substituting Eqn. (22 1 into 
Eqn. (pjl) obtains 



dG* 
dsj 



1 



i-nr=i«; 



n+l 



(28) 



which, upon inspection of Eqn. (24 1, is valid for j = n as 
well. 

Stability of the fixed point X* requires that the Jaco- 
bian of Eqns. ( [sP , 



-ri 

"2 -^2 



d' 



(29) 



be negative definite or, since the determinant is the prod- 
uct of the eigenvalues, that 



< (-l)"det(J) 

71 n 



fc=i 



1=1 



(30) 
(31) 

(32) 



fc=i 



1=1 



Since the decay rates are positive, Eqn. (32 1 says that 



the term inside the brackets in Eqn. (28 1 is positive for 
stable fixed points. 

For the networks in Fig. [l] where in the 1- and 2-cycles 
the reporter is attached by means of intermediates, the 



analog of Eqn. (28 1 is calculated similarly to be 



dG* 
dsj 



1=1 "/ 



N 
k=j+2 



(33) 



where = 4 is the number of genes, 1 < j < — 1 
for each of the 3 possible small molecule inputs, and n 
is the length of the cycle (1 < n < N — 1). Here 6 is 
the Heaviside function, for which we use the convention 
9{0) = 1. Its presence reduces the bracketed term to 1 
when the input Sj is outside the cycle, leaving only the 
contribution corresponding to the cascade from Sj to G, 
as must be the case. 



In Eqn. (33 1, the term outside the brackets represents 



the direct (i.e., the shortest) path from Sj to G and fixes 
the sign of dG*/dsj (since the term inside the brackets 
is positive at a stable fixed point). If the creation rates 
are monotonic (which is the usual model for transcrip- 
tional regulation, but may be violated in protein signal- 
ing due to competitive inhibition and other effects), this 
sign is unique and fixes the sign of AG* /Asj via Eqn. 
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Importantly, this says that the feedback in each of 
the topologies in Fig. [T|is irrelevant in determining the 
sign of AG* /Asj for a steady-state analysis. As an ex- 
ample, for the network in Fig. (inset), G* changes 
with increasing si according to a2a'^a'^, which, since Si 
inhibits the activation, is negative x positive x negative 
= positive, just as one would expect if the feedback was 
ignored. 



Direct functionality corresponds to specific 
orderings of output states 

Consider the case in which there are only two small 
molecule inputs. Si and S2, as in Fig. (inset). Since 
each input can be absent or present, 81,82 S 
there are four chemical input states c — 8182 G 

{ , — — ,++}. Direct functionality admits only 

two orderings of the four output states G*, and hence 
the functionality of the network is severely limited by its 
serial topology. To see this, note that for Fig. |2|i (inset) 
we have 

AG*/Asi >0 G;_ > G*__ and 

Gl^ > G*_^; (34) 



AG*/As2 >0 G% > G*__ and 

g;+>g;_. 



(35) 



These conditions permit only the following output order- 
ings, irrespective of biochemical parameters: 



Gl_ < G1+ < G;_ < G;+ or 

Gi_ < g;_ < G1+ < g;+. 



(36) 



These two orderings nevertheless allow the realization of 
a significant subset of all possible logical functions that 
one can build with two binary inputs, depending on the 
distinguishability of the four output states, as described 
in the next section. Quantifying the distinguishability 
demands careful treatment of the noise with a stochas- 
tic equivalent of our deterministic dynamical system, as 
described in the Appendix. 



NUMERICAL RESULTS 



We numerically solved the system in Eqns. (3][5l [with 
stochastic efi'ects given by Eqn. (6I1] with many param- 
eter settings for all 24 networks represented in Fig. [T] In 
addition to verifying the restriction to direct function- 
ality, we find that all networks can achieve all possible 
direct functions, suggesting that the networks are still 
quite versatile within the functional constraint. 

For all networks, we consider the case of two small 
molecule inputs Si and S2, as in the experimental setup 
of Guet et al. [T], and as shown for an example network 



Parameters 


Range 


decay rates, ri 


10""* - 10^^ 


Michaelis-Menten constants, Ki 


10° - 10^ 


basal expression levels, ao,i 


10"^ - 10"^ 


expression level ranges, at 


10° - 10^ 


cooperativities, hi 


10° - 10^ 


"on" input factors, 


10^ - 10^ 



TABLE I: Parameters and ranges from which each is ran- 
domly drawn, with 1 < i < 4 for the four genes, and 1 < j < 2 
for the two small molecule inputs. Ranges are representative 
of typical cell conditions (7| il8j . 



in Fig. (inset). We take sj to be a multiplicative 
factor by which the transcription factor concentration Xj 
is effectively scaled, i.e. 



(37) 



Then sJ = 1 for the "off" settings, and the sJ > 1 are 
free parameters for the "on" settings. 

We model the regulation using the familiar Hill form 
(which is monotonic and thus satisfies the direct func- 
tionality conditions) 



"(x) = Qo + Q ^;^^ (up-regulating), (38) 

a(x) = ao + a _ (down-regulating), (39) 

A" -I- X 



with basal and maximal expression levels Oq and ag + a 
respectively, Michaelis-Menten constants K, and coop- 
erativities h. For the 4-gene networks in Fig. [T] with 
only two small molecule inputs Si and S2, this gives 22 
parameters in total (cf. Table |l| . 

For a given parameter set, we numerically solve Eqns. 
([3]|5| (using Matlab's ode 15s) for each input state c G 
{ , — |-,H — ,++} to find mean steady-state concen- 
trations G*. We then solve Eqn. (61) to find fiuctua- 



tions around these means, giving probability distribu- 
tions P{G*\c) (cf. Appendix). The function is defined 
by the ranking of the conditional distributions P{G*\c). 
That is, if two distributions are distinguishable, then the 
one with the larger mean is ranked higher. We con- 
sider two distributions to be indistinguishable when their 
means are separated by less than the smaller of their stan- 
dard deviations (alternative definitions do not change our 
results qualitatively) , in which case they both take on the 
average of their two ranks. When there are only two dis- 
tinguishable output states, this rank-based classification 
reduces to that defining the familiar binary logical func- 
tions AND, OR, XOR etc. (see, for example, Fig.[2]3, (ii-v)). 
More generally, for one, two, three, and four distinguish- 
able responses, there are 75 total rankings (as listed on 
the horizontal axis of Fig. l2k). However, only 12 of these 
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satisfy the ordering constraints for each network analo- 



gous to those in Eqn. (36 1 and therefore correspond to 
direct functions (for the newtork in Fig. these 12 are 
shown in green on the horizontal axis). 

We ran 50,000 trials for each of the 24 networks, in 
which the parameters were randomly selected (using a 
distribution uniform in log-space) from the ranges in Ta- 
ble H] We found the steady-state reporter expression dis- 
tributions P{G*\c) and classified the responses by rank- 
ing. All 24 networks displayed only direct functions. 
However, every network was able to achieve all 12 of 
its direct functions with parameters selected via Table 
[l] meaning that the networks fully realized all the func- 
tionality allowed by the constraint. This suggests that 
the networks studied are both constrained and versatile, 
and that a cell may still use a serial network to perform 
multiple logical functions by varying biochemical param- 
eters, despite the restriction to direct functionality. Fig. 
[2] shows a histogram of functions and an example of each 
type of direct function for a representative network. 

We note that Guet et al. experimentally observed both 
direct and indirect functions [T . However, they explicitly 
call the indirect functions into question, citing several 
possible unanticipated effects including RNA polymerase 
read-through. We have not incorporated such effects into 
the current model. 



MULTIPLE FIXED POINTS 

For the 12 networks in which the overall sign of the 
feedback cycle is positive, there are parameter settings 
that support multiple stable fixed points. In this section 
we evaluate the extent to which the presence of multiple 
fixed points affects the constraint to direct functionality, 
and we find that violation of the constraint is possible 
but unlikely. 

While the function of a network has been defined in 
terms of P{G*\c), the linear noise approximation (cf. Ap- 
pendix) only gives us access to F( G* I c, X„ ) , the distribu- 
tion expanded around a particular fixed point X*^. The 
two are related by a weighted sum. 



P(G*|c)= ^^„P(G* icx;,). 



(40) 



where the probabilities tt^ of being near the mth fixed 
point will depend on the basins of attraction and cur- 
vatures near the fixed points. Numerical solution for 
P(G* |c) directly is possible in principle, although compu- 
tationally difficult. Whether the statistical steady state 
distribution is calculated numerically or is approximated 
as in this manuscript, if we continue to define the func- 
tion of the network by the ranking of the means of the 
P(G*|c), we have 



dG* 



d 

ds.j 



dG*G*P{G*\c) 



(41) 



Ranking 



(a) 



(i) Ranking: 2.5, 2.5, 2.5, 2.5 (li) Ranking: 1 , 3, 3, 3 (OR) (iii) Ranking: 2, 2, 2, i (AND) (iv) Ranking: 1 .5, 1 .5, 3.5, 3.5 (A) 




.5, 3.5, 1.5, 3.5 (B| (vi| Ranking: 1 , 2, 3.5, 3.5 




(x) Ranking: 1.5, 3, 1.5, 4 



{xi] Ranking: 1, 2, 3, 4 



Steady-State reporter concentration, G (proteins/cell) 
(b) 

FIG. 2: Direct functionality in a representative network with 
serial regulation (network shown in the inset in (a)), (a) His- 
togram of logical functions, as defined by the ranking of the 
output distributions P(G*|c) (cf. Numerical Results). Binary 
logic names are included after rankings when applicable, with 
'A' and 'B' corresponding to inputs Si and S2 respectively. 
Direct functions are labeled in green, indirect functions in 
red. Note that only direct functions are observed, and that 
all direct functions can be attained by the network, (b) An 
example of each direct function. Two distributions are consid- 
ered indistinguishable in rank when their means are separated 
by less than the smaller of their standard deviations. For ex- 
ample, in (ix), the distributions for the first and second states 

( and — h respectively) are indistinguishable, so they are 

tied in rank at 1.5 (the average of ranks 1 and 2). Note that 



all functions satisfy the ordering constraints in Eqn. ( 36 1 



Ed 



El dG*^ djTrn^* 
I T^rn— 1 ; — G, 



dSj 



dG*G*P{G*\c,X*J (42) 
(43) 



dsi 
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The expressions for the individual dG^^/dSj are given by 
Eq. (33 1, so the first term in Eqn. (43 1 exhibits direct 
functionahty. If the weights do not depend appre- 
ciably on the Sj, the second term will be small, and the 
restriction to direct functionality will be maintained. If, 
on the other hand, the weights do change appreciably (an 
obvious case might be the presence of a bifurcation at a 
particular value of Sj), then the second term may over- 
power the first enough to change the sign of AG* /As and 
violate the restriction to direct functionality. 

We investigate this effect in two ways. First, we show 
analytically that, in the case of a 1-cycle, crossing a bi- 
furcation does not violate direct functionality. Second, 
we subject all positive-feedback networks to a numerical 
test to estimate the dependence of the weights tt™ on 
the Sj. The results of both techniques suggest that the 
likelihood of a violation of direct functionality due to the 
presence of multiple fixed points is low. 



Bifurcations do not violate direct functionality (1-D) 
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Consider the case of a positive 1-cycle with a gene G 
immediately outside, as shown in Fig. [3^ (inset). For 
n = 1, Eqns. ( T0p2 1 become 



X* = ai(x*), 



(44) 
(45) 



where unnecessary subscripts are dropped and x~ X j s 
as in Eqn. (37). With a\ of the form in Eqn. ( [38| , there 
are at most two stable fixed points X\ and Xo, with 
X\ < X2 , as illustrated by an example in Fig. As 
shown in Fig. [3]d, bifurcations occur at si and S2 such 
that only X2 exists when s < Si, only X^ exists when 
s > S2, and Xi and X2 are found with (unknown) prob- 
abilities 7i'i(s) and 7r2(s) = 1 — •7ri(s) respectively when 
si < s < S2. These statements can be combined such 
that 

TTi{s) = 9{s - si)e{s2 - s)ni{s) + e{s - S2) (46) 

and 7r2(s) = 1— 7ri(s) define the probabilities of approach- 
ing Xi and X2 respectively for any s. Here 9 is the 
Heaviside function. 

As we go from an "off" value s" to an "on" value s+, 
let us assume that we hit both bifurcations, such that 
s~ < si < S2 < s+. To test for direct functionality, we 
investigate the sign of 



AG* _ 1 
As ^ As 



ds 



" m 

(48) 



A 

T1 + T2 



X 

(0 



(Kr/a)s 

(b) 

FIG. 3: (a) Solid: plot of regulation function ai — ai/r 
(refer to inset network), as defined in Eqns. (38 1 and (37 1, 
with parameters h = 2, ao = 0.03, a = 10, A" = 2, and r — 1. 
Dashed line shows ai = X such that dotted lines indicate 
locations of stable fixed points XI and X2 (take X| > Xl). 
(b) Stable fixed points X^ (solid) and unstable fixed points 
(dashed) as a function of s, with ao/a — 0.003 as in (a). 
Dotted lines indicate locations of bifurcation points si and S2 
such that only X| exists when s < si, only X* exists when 
s > S2, and both Xi and exist for si < s < S2. 



obtained using Eqns. ^ and The first term Ti 

depends on 



dGl 
ds 



a2 



1 - a[ 



(49) 



(from Eqn. (28 1; a' = da/dX and a = da/ds as before, 
both evaluated at the mth fixed point), which, as pre- 
viously discussed, is always of the sign of 012, consistent 
with direct functionality. 
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The second term T2 can be written 



1 
1 



dni dTT2 , 

as as 



diTi 
ds 



(g; - Gi) ds, 



(50) 
(51) 



and since 



diTi 
ds 



^(s2 - s)tti{s)S{s - Si) 

+ [l-0(s-si)niis)]6is~S2) 



+9{s - si)9is2 ~ s) 



dfti 
ds 



(52) 



(where 6 is the Dirac delta function) we have 

T2 = ^{-[^i(G;~GD],^-h(G;-GD], 
dni 



ds 



iG*2~Gl)ds 



(53) 



The first two terms in Eqn. (53 1 represent the contribu- 



tions from crossing the bifurcations at si and S2 respec- 



tively. Using Eqn. (45 1 we may write them as 



To = 




Aa2 



{G*2~Gl)ds 



Ax* 



(54) 



where Aa2 = a2(X2) - o:2iXi) and Ax* = X2 - Xi = 
{X2 — Xi)/s > 0. Since a2 is monotonic in X, 
— Aq:2/Ax* at fixed s is of the same sign as —a2, which 
is of the same sign as 6(2 since s effectively reduces X 
(Eqn. (37)). Therefore the contributions to AG* /As 



from crossing the bifurcations do not violate direct func- 
tionality. A violation, at least in the case of a 1-cycle, 
can only come from variations in the probabilities iTm 
within the region si < s < S2, as described by the last 



term in Eqn. (54). Next we describe a numerical test 



that suggests such violations are rare. 



Numerics suggest violations from multiple fixed 
points are rare 

For each of the 12 positive-feedback networks, we nu- 
merically found the steady state of the dynamical system 
with randomly sampled parameters as before (cf. Numer- 
ical Results). However now for each parameter set we 
solved the system many times with randomly selected ini- 
tial conditions. When multiple fixed points were found, 
the fraction of trials approaching the mth fixed point 
was used for the weight iTm- This assumes the iTm are 
determined only by the basins of attraction of each fixed 



point, and by the distribution of the initial conditions. 
However, different distributions of initial conditions do 
not result in qualitative different results. 

For each network, 2,000 parameter sets were selected 
(uniform randomly in log-space), at which the system 
was solved 100 times with initial protein counts selected 
uniform randomly from to 1, 000 proteins per cell. Over 
all positive-feedback networks, 37% of the parameter sets 
supported multiple fixed points for at least one of the set- 
tings of Si and S2. However only 0.46% of parameter sets 
produced violations of direct functionality. Moreover this 
number is likely an overestimate, as no distinguishabil- 
ity criterion was imposed as was done in the single-fixed 
point case (cf. Numerical Results). It is likely that this 
fraction would remain low if the estimation of the tt^ was 
refined to incorporate the curvatures of the fixed points, 
or if alternative distributions were used for the sampling. 



ALL NETWORKS WITH SERIAL REGULATION 
EXHIBIT ONLY DIRECT FUNCTIONALITY 

In this section, we extend our analytic constraint as 
derived in the context of the system studied experimen- 
tally by Guet et al. [1] to show that any network with 
only serial regulation — each node having or 1 parent — 
exhibits only direct functionality, i.e. any target node 
changes with any input Sj according to the direct path 
between them. 

We first consider a connected directed graph in which 
every node has in-degree 1, called a contrafunctional 
graph |12j . One can show that a contrafunctional graph 
has exactly one cycle, each of whose nodes is the root of 
a tree if the cycle edges are ignored [12]. Now consider 
changing one node's in-degree to 0, or equivalently, re- 
moving an edge. If the edge is in the cycle, the graph 
remains connected and becomes a tree. If the edge is not 
in the cycle, the graph is cut into two components: a 
contrafunctional graph and a tree. 

A tree exhibits only direct functionality since there is 
at most one path from an input Sj to a gene X^, which 
is therefore the direct path. 

In a contrafunctional graph, we first consider the case 
where the target node X^ is inside the cycle. Only inputs 
Sj that are inside the cycle can affect Xj because the rest 
of the graph consists of trees that all point away from the 
cycle. Since we can start labeling nodes at any point in 
the cycle, we may take i < j without loss of generality. 
Then, using the chain rule. 



dxi 

dsj 



dx; 



'■{j modn) + l "f' '^{j modn) + l 



n 



k=i ' 



n 



(55) 



l=i "(;mod n) + l 
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n 



fc=i ' 



n 



l=i '^(imodn) + l 



(56) 



where the second step follows from Eqn. (22) 



We next consider the case where the target node is 
outside the cycle. An input Sj can only affect the node 
if it is either in the cycle or above the node in its tree. 
The portion of the path in the tree will exhibit direct 
functionality. Therefore in looking for possible indirect 
functionality we may, without loss of generality, take the 
node to be immediately outside the cycle, as we did for 
G in the previous section. dG* /dsj is then given by Eqn. 



(281 



In both Eqns. (56 1 and (p8|, the term outside the 



brackets represents the direct path from Sj to the tar- 
get node, and the term inside the brackets is positive for 
stable fixed points. Therefore, a contrafunctional graph 
exhibits only direct functionality. Since each connected 
component of a network in which every node has in- 
degree or 1 is either a contrafunctional graph or a tree, 
such networks exhibit only direct functionality. Thus, in 
general, the possible logical functions of topologies with 
at most one regulator per node are severely constrained. 



APPENDIX: THE STOCHASTIC MODEL 



The dynamical system in Eqns. (3|5) provides a deter- 
ministic description of mean expression levels but fails to 
capture fluctuations around these means. A full stochas- 
tic description is given by the chemical master equation. 
For N species participating in R elementary reactions in 
a system with volume 17, the master equation reads 



dP{n,t) 
dt 



R 

■E 



N 



1 /,(X,f7)P(n,t), 



(57) 

where P(n, t) is the probability of having the copy num- 
ber vector n = fiX = ri(Ai, . . . , Ajv) at time t, Zij is the 
N X R stochiometric matrix, E~^'^ is the step operator 
which acts by removing Zij molecules from n^, and fj is 
the rate for reaction j. The fj are the ctj and rjXj of 
Eqns. (3|5l in the macroscopic limit. 

As in previous work 0, we employ the much-used lin- 
ear noise approximation [T31 HH HSl US] to make Eqn. 
(57 1 tractable by expanding in orders of fi"^/^. Intro- 
ducing ^ such that rii = fiA^ -I- fi^^^^i and treating ^ as 
continuous, the first two terms in the expansion yield the 
macroscopic rate equations (e.g. Eqns. (3]|5l in our case) 
and the linear Fokker-Plank equation, respectively: 



N 

E 

i=l 



dt d^. 



N R 

EE^> 

i=l 3 = 1 



./.(X)^, (58) 



dP{^,t) 
dt 



(59) 



__^z,j{df,/dXk) 



where Jjfc = 

trix (e.g. Eqn. (^291) and Di^ — X]j=i 
diffusion-like mamx. The steady-state solution to Eqn 



is the Jacobian ma- 
R 



^ zjij Zkj fj(X.) is a 



( 59 1 is the multivariate Gaussian 



P(0 = [(2^)^detS] '/%xp 



where the covariance matrix ^ satisfies 



JS + S 



D = 0. 



(60) 



(61) 



We solve for S using standard matrix Lyapunov equa- 
tion solvers (e.g., Matlab's lyap). Thus fluctuations are 
captured to leading order by Gaussian distributions with 
means given by the macroscopic equation and vari- 
ances given by the diagonal entries of S. For example, 
Gaussian distributions P{G*\c) are shown in Fig. [2jD for 
the steady-state concentration of the reporter gene G un- 
der chemical input states c. In [6] we have compared the 
distributions P{G*\c) obtained using the linear noise ap- 
proximation to those obtained via direct stochastic simu- 
lations [17 and found the results almost indistinguishable 
for molecular copy number above 10-20. 
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